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I. INTRODUCTION 

\ The nonlinear Schrodinger equation (NLS) is one of the universal models of mathematical physics arising in many 

■ different physical scenarios [ll, Q- In recent years there has been an enormous interest on the study of NLS equations 
with nonlinear coefficients depending either on space, time or both. The motivation comes from the applications of 
the model to the fields of Bose-Einstein condensates (BECs) and nonlinear optics. In the first field of application, the 
so-called Feschbach resonance management allows for a precise control of the atomic interactions responsible for the 
strength of nonlinearities. This has been experimentally implemented leading to time dependent nonlinear coefficients 
and would allow for the generation of space or space- time dependent nonlinearities. 

Many different types of nonlinear waves have been studied theoretically, experimentally or numerically in the 
context of BECs (see e.g. the reviews [sl, However the consideration of inhomogeneous nonlinearities has led to 
^sq \ the prediction of many remarkable nonlinear phenomena in the last few years, either for time dependent (5l-[l5| or 

. space dependent nonlinear coefficients p!6l-[37|. 
CN| ■ Although the model equations of interest are not integrable in general, there have been efforts to construct families 
of exact solutions of the NLS with spatially inhomogeneous coefficients. Specifically, grou p-theoretical methods based 

■ on Lie symmetries have been used to obtain solutions of NLS type equations in Refs. |24 l25l, [3ll, [36l. IssMioj . However, 
\ up to now the applications of the technique have been restricted to scalar systems. The purpose of this paper is to 

extend the method and construct exact solutions for multicomponent systems. These systems are of an enormous 
physical interest because of the extra degrees of freedom provided by the existence of two independent subsystems, 
. e.g. in term of application to multicomponent BEC the model describes two hyperfine levels of an atomic BEC. Exact 
^\ ■ soliton solutions for such systems in the case of homogeneous in time and space nonlinearities have been reported in 
many papers (see e.g. [41r.-47] and references therein). Very recently in [48| the dynamics of the multicomponent BEC 
with spatially modulated nonlinearity was treated by means of the the variational approach. Here we will go beyond 
previous studies and use Lie symmetry method to construct exact bright-bright, dark-dark and dark-bright solutions 
of the two-component system with spatially inhomogeneous nonlinearities and investigate their linear and dynamics 
stability. 

The paper is organized as follows. First, in Section [III we introduce the general theory of the Lie symmetry analysis 
for our model problem: a system of two coupled nonlinear Schrodinger equations with spatially inhomogeneous 
nonlinearities. In Section Hill we use the method to construct explicit solutions for systems without external potentials 
{Vj{x) = 0, j = 1, 2) and study the linear and dynamical stability of the obtained solutions. In Section [IVl we present 
some families of exact solutions for inhomogeneous coupled nonlinear Schrodinger equations with external potentials 
quadratic in space and propose their stability analysis. Next, in Section |Vl we discuss the construction of dark-bright 
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soliton solutions and study their stability. Finally, in Section IVI) we summarize our conclusions. The theory of the 
linear stability analysis for the system is presented in Appendix \K\ 

II. GENERAL THEORY 
A. Model to be studied 

In this paper we will construct and study the stability of solutions of the following coupled spatially inhomogeneous 
nonhnear Schrodinger equations (CINLS) 

= -^ + Vi{x)i^, + {gn{xMi\^+gu{x)\H^)i^i, (la) 

= -^ + V2{x)i^2 + {g2l{x)\i^l\^+g22{x)\H^)i^2, (lb) 

for known real potential functions Vj{x), and spatial modulations of the nonlinearities described by the real functions 
Qjkix)^ j, /c = 1, 2. In applications to BEC, ijji and ^1)2 are complex wave functions describing two hyperfine levels (the 
extension to different atomic species is simple, corresponding to the addition of different masses) and describe the 
dynamics of coupled quasi-one dimensional BECs in the mean field limit with spatially inhomogeneous interactions. 

In the following we will be interested in solitary wave solutions of Eqs. ([T]) of the form t) = Uj{x)e~^^^'^ ^ where 
Uj{x) are real functions (j = 1,2). Upon substitution in Eqs. ^ we get system of stationary equations 

^lUi = ^Vi{x)ui ^ {gii{x)ul^ gi2{x)ul) ui, (2a) 

A2^2 = + V2{x)u2 + {g2i{x)ul + g22{x)ul) U2. (2b) 

The construction of families of exact solutions of Eqs. ([2]) is the main goal of this paper. 

B. Lie symmetry method for Eqs. ([2|) 

In this section we will extend the method of Ref. [2^ to the CINLS given by Eqs. (j2j). The methodology is similar 
although the calculations are more complicated, but we will show that it is still possible to extend the method for 
this more complicated situation. 

It is said that a second-order differential system 

Nj{x^Uj^Uj_^xi'^j,xx) — O5 (3) 

possesses a Lie point symmetry [49| of the form 

d d d 
M = ^{x,ui,U2)— ^r]i{x,ui,U2)- hr]2{x,ui,U2)^—, (4) 

ox OUi OU2 

if the action of the second extension of M, i.e. M^^^ on Nj^ gives zero whenever (|3]) is satisfied, i.e. 

r 



^[X,Ui,U2) — ^Vl{XiUi,U2)- 

OX OUi 



r]2{x,ui,U2)- hr^i {x,ui,U2,uix)- h r/2 (^, ^1, ^2, ^2,^)- 



OU2 OUi^x OU2,x 

Nj{x,Uj,Uj^x,Uj^xx) = (5) 



d 

-f T^f ^ (X, Ui,U2, Ui^xx)^ h r?2^^ (X, Ui,U2, U2^xx) 



I ^ \ • ri y^i ^z,xxj ^ 

OUi^xx OU2,xx 



where Uj^x = duj/dx and Uj^xx = d^Uj/dx'^^ with j = 1, 2. 
In our problem, Nj{x^Uj ^uj^xx)^ j = 1^2 are given by 



Ni{x,ui,U2,ui^xx) = ui^xx - h{x,ui,U2), (6a) 

N2{x,Ui,U2,U2^xx) = U2,xx - f2{x,Ui,U2), (6b) 
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where 





0, 


S,UiU2 ~ 


0, 


£,U2U2 ~ 


0, 


Vl,U2U2 


= 0, 


Vl,UiU2 


~ ^U2,x 5 


Vl,U2X - 


~ S,U2fl^ 




'^^XUi 1 




^xx 


r]i,xx + 





h{x,ui,U2) = \Vi{x) - gii{x)ul^ gi2{x)ul]ui, (7a) 
f2{x,ui,U2) = [V2{x) - X2^ g2i{x)ul^ g22{x)ul]u2, (7b) 

and the action of the operator M^^^ on A^i leads to 

(8a) 
(8b) 
(8c) 
(8d) 
(8e) 
(8f) 
(8g) 

a^nJi-Mn, =0, (8h) 
2^x) fl - ifi,x - Vifi,ui - V2fi,u2 + f2Vi,u2 = 0- 

In a similar way, we get equations for the action of M^^^ on N2 (we omit the details). Solving the previous equations, 
we find that the only Lie point symmetries of Eq. (j2j) are of the form 

Odd 
^ c(x)— + m{x)ui- h m{x)u2- — , (9) 

ox OUi OU2 

where 

9jk{x) = gjk/c{xf, (10a) 
m{x) =c'{x)/2, (10b) 

-c{x)Vi{x) + - 2c' {x) (Fi(x) - Ai) = 0, (10c) 

V2{x) = Vi{x) + K/c\x) + A2 - Ai. (lOd) 

Here gjk =constants and K an arbitrary real number. 

It is known [50|, that the invariance of the energy is associated to the translational invariance. The generator of 
such a transformation is of the form M = d/dx. To use this fact, we define the transformation 

X = h{x)^ Ui = v{x)ui^ U2 = v{x)u2^ (11) 

where h{x) and v{x) will be determined by requiring that a conservation law of energy type M = d/dX exists in the 
canonical variables. Using Eqs. ([9]) and ([TT]) . we get 

h{x)= f -^ds, v{x) =c(x)-^^^. (12) 
Jo c{s) 

Rewriting Eqs. (|2]) in terms of Uj = c~^^'^{x)uj , j = 1, 2 and X = l/c{s)ds we get 

- ^ + {hiU! + ~gi2Ui) Ui = E,Ui, (13a) 



d'U2 

dX^ 



-^{g2iUl^~g22Ul)U2 = E2U2, (13b) 



where 



El = {Xi-Vi{x))c{x)^ -lc\xf ^lc{xy{x), (14a) 

\2 



E2 = {\2-V2{x))c{xf-lc'ixf+'^cix)c"ix), (14b) 

are constants. Thus, in the new variables we obtain nothing else but system of NLS equations without external 
potential and with spatially homogeneous nonlinearities. Of course not any choice of Vj{x) and gjk{x)^ j^k = 1,2, 
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leads to the existence of a Lie symmetry or an appropriate canonical transformation since the function c{x) must be 
sign definite for /7j, j = 1,2 and X to be properly defined. In the specific case when all nonlinear coefficients are 
equal, i.e. gjk = g for all j, /c = 1,2, Eos. ([T3j) correspond to the time independent Manakov system [4l|, that has 
been studied in detail (see, for instance [42-44]). 

We can thus use the known solutions of Eqs. ()14b|) to construct solutions of Eqs. ([2]) in the same spirit as in the 
scalar case (see e.g. Ref. [223). Moreover, when both Ei = E2 = E and the restriction 

(^11 - ^2l)(^22 - 912) > 0, (15) 

are satisfied, Eqs. (fT3|) reduce to the scalar NLS equation 

-^+a^'=E^, (16) 

where we have introduced the new field <l>, related with the original fields through the relation 

Uj=Oj^. (17) 



Here Oi = ^1^22 -^i2|/Q, and O2 = \/\hi^^^2i\/Q , being Q = \gii922 - 9i292i\ and a = sgn[{giig22 - 912921) {9 jj - 
912)] which due to the condition (fT5]) does not depend on j. In that situation we can use the many known solutions 
of the scalar equation (p!6]) to construct solutions of the vector problem (|2]). Equation (p!6|) has many known exact 
solutions. Remarkable examples, to be used later in this paper, are 

^ = /J. tanh (^-^^^ ' ^ = V, ^ = 1' (18a) 

^ = V2fi ] E = -^\ a = -l, (18b) 
cosh(/iA ) 

^ = My2(r^£|^|^, E = ^\l-2k'), a = -l, (18c) 

<!> = V2/idn(/iX,A;), E = ii^{k^-2), a = -1. (18d) 



III. SYSTEMS WITHOUT EXTERNAL POTENTIAL {Vi{x) = V2{x) = 0) 

A. General considerations 



As a first application of our ideas let us choose Vi{x) = V2{x) = 0. By taking K = and A2 = Ai = A in (|lQd|) we 
find from Eq. (|10cp the equation c^^\x) + 4Ac'(x) = 0, whose general solutions are given by 

c{x) = Ci sin cjOx -\- C2 cos ux -\- Cs (A > 0), (19a) 
c{x) = Cie'^^ + Cse-^^ + Cs (A < 0), (19b) 

where uj = 2Y^^lAj. It should be emphasized that the auxiliary function c{x) is a periodic function for A > which 
corresponds to the continue spectrum of the linear eigenvalue problem cPm/dx'^ + 4 Am = and for A < (the 
region where the linear eigenvalue problem does not have periodic solutions) the function c{x) becomes exponentially 
localized. Therefore, in the following subsections we will distinguish two types of solutions of Eqs. ([1]): corresponding 
to either periodic or localized nonlinearities. 



B. Dark-dark soliton solutions 



Let us first consider the case A > 0. Then, taking Ci = 0, C2 = a, C3 = 1, Eq. ()19ap leads to a periodic dependence 
of the coefficient, e.g. c{x) = 1 + acosux, and according to Eq. (jlQap the nonhnear coefficients are given by 

9jk{x) = 9 jk {I ^ a cos ujxy^ , j,k = 1,2. (20) 

Although the form of nonlinear coefficients seems complicated, it is easy to see that for small a this nonlinearity is 
approximately harmonic in space gjki^) — 9jk (1 — Sacoscjx), a <C 1, j, /c = 1, 2, and thus at least in that limit is of 
direct physical interest. 



FIG. 1: Dynamics of the dark-dark solution ()2ip under a perturbation of the form ([23]). Shown are (a) the initial density 
profiles of the first (black soUd) and second (red dashed) component. (b,c) Density plot of the time evolution of (b) |i^i(x,t)|^ 
(c) \u2{x,t)f. Parameter values are g — 0.1, = 1, ^22 = 0.5, a — 0.3, cj = 1, e = 0.001, q — 0.01. 

We can construct our canonical transformation by using Eqs. (fTT]) and obtain 

tan(-Vl-«2x(x)J = tan—. (21) 

Using any solution of Eq. (p!3]) with E = Ei = E2 = (l — o?^ > 0, where A = a;^/4, this transformation provides 
solutions of Eq. (|2j) with g{x) given by Eq. ([20]) . For instance, taking the nonlinear coefficients to satisfy Eq. ([T5]) 
and cr = 1 we can construct dark-dark soliton solutions of Eq. ([2]) of the form 

Uj{x) = Oj (1 + a cos(cjx))"'"^^ A/2/i tanh (/iX(x)) , (22) 

where /i = = ^uj\^1 — o? and X(x) is defined by Eq. ([2T]) . We want to emphasize that this is only a simple 
example of the many possible solutions that can be constructed in such a way. 

C. Stability of the dark-dark solution ([22]) 

Let us consider the dynamical stability of the dark-dark soliton solution ([22]) . Here and in all of our calculations to 
be presented in this paper in order to check dynamical stability we perturb the initial exact solution according to 

Uj(x) = ^ijo[l + ecos(gx ± 1)], (23) 

where and " correspond to the first and second components, respectively, e will be taken as a small number 
and q is the perturbation wave number. 

For the analysis we will fix the nonlinearity for the first component = 1 and change the nonlinearity for the 
second component in the interval < ^22 < 1 and the intercomponent interaction ^12 = ^21 = ^ in the interval 
— y/g22 < g < g22- Our results show that for small a <C 1 the dark-dark solution is dynamically stable in the full 
range < ^ < ^22 (see Fig. [1]). However, in the interval —^fg^/i < g < the dark-dark solutions are unstable (see 
Fig. [2]). 

Also we have checked that by increasing the amplitude a of the nonlinear inhomogeneity we found the instabilities 
take longer times to develop. As one can see in Fig. [3] for a = 0.9 the instability occurs near t ~ 250 while in the 
previous example for a = 0.3 shown in Fig. [2] the solution breaks around t ^ 90. 

D. Bright-bright soliton solutions 



Let us now consider the case A < defined by Eq. ([19bp . In order to simplify the analysis we restrict ourselves 
to two particular choices of the constants cj, Ci, C2 and C3. The first one is = 1, Ci = C2 = 1/2 and C3 = in 
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FIG. 3: Same as in Figfflbut for a = 0.9. 



Eq. (jl9b[) . Thus c{x) = coshx and A = —1/4, E = Ei = E2 = 1/4, being cosX{x) = — tanhx, and the nonhnear 
coefficients are 

gjk = — % — (24) 
cosh (x) 

Then < X < tt, and to meet the boundary conditions Uj{±oo) = 0, one has to impose Uj{0) = Uj{7r) = 0. This 
means that the original infinite domain in Eqs. (j2j) is mapped into a bounded domain for Eqs. ([13]). By choosing 
the nonhnear coefficients to satisfy Eq. ([15]) and <j = — 1 we can use the solution (|18cp to construct vector soliton 
solutions of the form 

f/,m = .,..^2(TT^£g||, ,25) 



where /i = 

The functions Ui{X) and U2{X) satisfy Ui{0) = U2{0) = and in order to meet Ui{7t) = U2{7t) = 0, the condition 

/iTT = 2nK(k) where K(k) is the elliptic integral K(k) = (l — k'^ s\v? Lp) d^p must hold. Thus, to satisfy the 
boundary conditions, k must be chosen to solve the following equation 



AnK (k) Vl - 2A:2 = TT for n = l,2,3,. 



(26) 
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FIG. 4: Dynamics of a bright-bright sohton solution ()27p with n = 1 under a perturbation of the form (|23|) . Shownare (a) the 
initial density profiles of the first (black solid) and second (red dashed) component. (b,c) Density plots of the time evolution 
of (b) \ui(x,t)\'^ (c) \u2(x,t)\'^. Parameter values are g — 0.2, gn — —1, ^22 = —0.5, e — 0.001, and q = 0.2. 




FIG. 5: Same as in FigHbut for g = 0.23. 

It can be shown that for every integer number n this algebraic equation has only a solution kn which means that 
there are an infinite number of solutions of Eqs. ([2]) of the form given by Eqs. (|25]) . Moreover, each of those solutions 
has exactly n — 1 zeroes. Then, using (fTTj) we find a solution of Eqs. ([2]) of the form 

Uj{x)=c{x)'^%{X{x)). (27) 

E. Stability of the bright-bright solution ([27|) 

Next we investigate the stability of the bright-bright solutions given by Eq. ([27|) by fixing gu = —1 and choosing 
diflFerent combinations of (^22,^) (as in the previous case we consider gi2 = ^21 = g)- First we start with the case 
n = 1 in Eq. ([26|) . In this case for each value of ^22 we find that there exists a critical value of ^ = gcr such that for 
g < gcr the vector solution is stable and for g > gcr the solution becomes modulationaly unstable. As examples Figs. 
[3] and [5] show that the critical value gcr lies between ^ = 0.2 and ^ = 0.23. It should be stressed here that in both 
regions of existence of the solutions, i.e. g < gcr and g > gcr^ they are linearly stable, i.e. there are no negative or 
complex eigenvalues l^j, {j = 1,2). However dynamically, the solutions for g > gcr suffer a phase "separation" where 
one of the component extrudes the other component, a phenomenon also observed in related systems Ref. (sij . 

As to solutions with n = 2 and 3, the linear stability analysis shows that all high-order solutions are linearly unstable. 
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• (c) 



Q.5 



(e) 



FIG. 6: Dynamics of a n = 2 bright-bright solution given by Eq. ([27|) under a perturbation of the form ([23|) . Shown are (a) the 
initial density profiles of the first (black solid) and second (red dashed) components. (b,d) Density plot of the time evolution of 
(b) |i^i(x,t)|^ (d) \u2{x,t)\'^. (c,e) Part of the spectrum of the eigenvalues ^^(Qr) for the first (c) and second (e) components, 
calculated from Eq. ()A5|) . The red dots on the spectrum show the existence of negative and complex eigenvalues. Parameter 
values are g — —0.1, gn — —1, ^22 — —0.5, e — 0.001, q — 0.2. 






(c) 










• (e) 




• 
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FIG. 7: Same as in Fig. [6] but for g = 0.1. 



In the case n = 2, the linearized equations (see Appendix [Aj) have one pure negative eigenvalue Qr = Re(r^) < for 
both components and also two complex conjugated eigenvalues with non-zero imaginary part Qi = lm{Q) 7^ for the 
first component only (see panels (c), (e) in Figs l6] and [7|). 

The numerical integration of Eqs. ([1]) with perturbed initial bright-bright solution with n = 2 confirms the results 
of the linear stability analysis. As we can see in panels (b) and (d) of Figs. [6]and[71 the dynamics becomes unstable for 
small times and consists of oscillatory behavior at the initial stages of instability. When the instability develops fully 
the initial profile of the solution (having two peaks) degenerates into a one-peak stationary profile corresponding to 
bright-bright solution with n = 1 (Figl6]). In the case of ^ > due to the instability, the second component splits into 
two wavepackets which move outwards the center while the first component transforms into one-peak bright soliton 
solution (see Figd), that is a stable mode of the system. 

In the case of n = 3 one more pare of complex conjugated eigenvalues appears in comparison with the case n = 2 
what corresponds to the more pronounced instability which starts to develop earlier then in previous case and in the 
case ^ < consists of transformation of the three-peaks profile of the solution into one-peak oscillating solution (see 
Figsl8]and[9]) while in the case ^ > the perturbed bright-bright solution disappears after strong radiation emission- 



IV. SYSTEMS WITH QUADRATIC POTENTIALS 
A. Construction of families of solutions 



In this section, we construct explicit solutions of Eqs. (|2]) when either (or both) of the external potentials Vi or V2 
are different from zero. We will consider the case when Ai = A2 = A and Ei = E2 = 0. By choosing c{x) = e~^^ and 
substituting into Eqs. (fT4|) and (|14b[) we get quadratic trapping potentials of the form Vi{x) = V2{x) = X^x^. These 
Gaussian nonlinearities can be generated by controlling the Feshbach resonances optically using a Gaussian beam (see 
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FIG. 8: Dynamics of a bright-bright solution given by Eq. ([27|) with n = 3 under a perturbation of the form ([23|) . Shown 
are (a) the initial density profiles of the first (black solid) and second (red dashed) component. (b,d) Density plot of the time 
evolution of (b) \ui{x,t)\'^ (d) \u2(x,t)\'^. (c,e) Part of the spectrum of the eigenvalues Qi(Qr) for the first (c) and second (e) 
components, calculated from Eq. (|A5p . The red dots on the spectrum show the existence of negative and complex eigenvalues. 
Parameter values are g — —0.4, gn — —1, and ^22 = —0.5, e = 0.001, q = 0.2. 
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FIG. 9: Same as in Fig. |H]but forn = 3 and g = 0.4. 



e.g. [S2l)- Here we also take K = 0, thus 

9jk{x) = gjk exp (SAx^) , j,k = 1, 2. 



(28) 



J erf \/—Xx. This leads us to a restriction 



Our canonical transformation is given by X{x) = ds exp (As^ 

on the sign of the chemical potential A which must be less then zero, A < 0. In this case Eq. becomes 

(PUi 



' dX^ 



- giiU! + gi2UiUi = 0, 



(29a) 
(29b) 



? < X < i-v/— J, and hence, we can again construct many local- 



Note that the range of X is again finite since — a— — 2V a 
ized solutions to Eq. ^ starting from solutions of Eqs. (|29)) which satisfy the boundary conditions {7j(±iy^— ^) = 0, 
J = 1,2. 

Using the condition (fT5|l we can reduce Eqs. ()29|l to 



dX^ 



cr$"^ = 0, 



(30) 



where we introduce new field <l> which is related with original fields through the relation (fT7|) and a = sgn[(^ii^22 
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9i292i){9jj — ^12)]- For a = —1 this equation has two solutions 

^^^\X) = /icn(/iX,/c*), (31a) 

^ ^ V2dn(/iX,A:*)' ^ ^ 

with /c* = l/V^, which solve the system (j29l) and that <I>'^^)(X) and ^^'^\X) vanish when /iX = (2n + l)i^(/c*) and 
fiX = 2nK (A:*) correspondingly. Thus we come to an infinite number of solutions of system (|29|) under zero boundary 
conditions on the new finite interval, which correspond to different values of /i. 
Finally, localized solutions of the NLS system (|2]), are given by 

Ujn{x) = 2n^^K{K)eje''^^^ cii{xn{x),K), n = l,3,... (32) 

and 



with 



Xn{x) = nK (k^) erf (34) 



It can be shown by simple asymptotic analysis that the last factors in Eqs. (|32]) and (|33|) tend to zero as x ^ ±00 
faster than exp(— x^/2) and that these are indeed localized solutions of our problem. 



B. Dynamics and stability of the solutions 



Let us now check the stability of the bright-bright soliton solutions in the presence of the parabolic potential. First 
we consider solutions given by Eq. (j32j) with n = 1. As in the case of the bright-bright solution without external 
potential from Section UlI DI here solution for n = 1 is also linearly stable in all the domain of existence of the solution. 
However, by checking the dynamical stability under small perturbation of the initial exact solution we found that for 
fixed values of intra-species interactions gu = — 1 and ^22 = —0.5 there is the domain of stability ^Jg\\g22 < ^ ^ (see 
an example in Fig. [TQ|) and also that there is a domain where the solution is unstable ^ > (an example in Fig. [TT]) . 
As it is clear from Fig. [TlT b). (c) there is a phase separation effect related to the immiscibility of the two-component 
mixture. In this case the fist component is shifted to negative x values while the second component is displaced to 
positive X values; however this separation rather small and components continue to have a significant overlapping. 

Now let us consider stability of the bright-bright solution for n = 2. In contrast to the case without linear potential 
studied previously, here the solution is linearly stable in all of its domain of existence. From the dynamical stability 
analysis we find that, in the domain ^Jgxxgri < ^ ^ 0, the solution is dynamically stable while for ^ > as in the 
previous case of the n = 1 solution becomes unstable due to the immiscibility condition between the two components. 
Examples are shown in Fig. [12] that has a stable dynamics and in Fig. [13] where unstable behaviour of the bright- 
bright solution for n = 2 is observed. In contrast to the case of bright-bright solution without external potential in 
the case at hand the second component due to refiection from the parabolic potential starts to oscillate. 

Finally, when moving to solitons with n = 3 we observe that in the full domain of existence of the bright-bright 
solution, it is linearly as well as dynamically unstable (see Figs. [14] and [T5]) . 



V. DARK-BRIGHT SOLITON SOLUTIONS 



A. Dark-bright soliton solutions 



Up to now we have focused our attention in constructing vector extensions of the solutions known for the scalar case 
with spatially dependent nonlinear coefficients and studying their stability. In this section we will consider nontrivial 
solutions where the functional form of both components is completely different. 

Dark-bright soliton solutions have been extensively studied in the context of BEC after the pioneering work of 
Busch and Anglin_j45|, where the coefficients gij = 1, j, /c = 1,2, what corresponds to the well-known integrable 
Manakov system |41| . 
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FIG. 10: Dynamics of a n = 1 bright-bright sohton solution ()32p in a parabohc potential Vj{x) = A^x^, (j = 1,2) under a 
perturbation of the form (|23p . Shown are (a) the initial density profiles of the first (black solid) and second (red dashed) 
component. (b,c) Density plot of the time evolution of (b) \ui(x,t)\'^ (c) \u2{x,t)\'^. Parameter values are g — —0.1, gn = — 1, 
^22 = -0.5, A = -1, e = 0.001, q = 0.2. 




FIG. 11: Same as in Fig. [TO]but for g = 0.4. 



As in the dark-dark case to constract dark-bright solution it is netural to consider Vi (x) = for the dark component, 
Ai = A2 = A = a;^/4 > and correspondingly c{x) = 1 + acosux. From Eq. (p!3|) we get the relation 

El = A(l-a^), (35a) 
E2 = Ei-V2{x)c{x)\ (35b) 

and taking the potential V2{x) of the form 

V2{x) =i^(l + acos(cJx))-^ (36) 

from (j35aj) we get that Ei 7^ ^2, and 



K = Ei-E2. 



(37) 
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FIG. 12: Dynamics of a n = 2 bright-bright sohton solution (|33)) in a parabohc potential Vj{x) = A^x^, {j = 1,2) under 
a perturbation of the form ([23]) . Shownare (a) the initial density profiles of the first (black solid) and second (red dashed) 
component. (b,c) Density plot of the time evolution of (b) \ui{x,t)f (c) \u2{x,t)f. Parameter values are g = —0.1, gn = — 1, 
^22 -0.5, A = -1, e = 0.001, q = 0.2. 
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FIG. 13: Same as in Fig. [IO]but for g = 0.2. 
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FIG. 14: Dynamics of a n = 3 bright-bright soliton solution (|32p in a parabolic potential Vj{x) — A^x^, (j = 1,2) under 
a perturbation of the form ([23|) . Shownare (a) the initial density profiles of the first (black solid) and second (red dashed) 
component. (b,c) Density plot of the time evolution of (b) |i^i(x,t)|^ (c) \u2{x,t)\'^. Parameter values are g — —0.4, gn = — 1, 
g22 = -0.5, A = -1, e = 0.001, q = 0.2. 
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FIG. 15: Same as in Fig[Tl]but for g = 0.4. 
Therefore, we can choose as solutions of Eqs. ([T3|) the fohowing dark-bright sohton solutions: 

Ui{X) -- 



U2{X) 



tanh(/iX), 

911 

I El {§21 -gii) 
gii {922 -912) 



sech(/iX), 



where 



' El {gng22 -912921) 
2^11 {922-912) 



E, = ^^Ei-f,'. 
911 



(38a) 
(38b) 

(39a) 
(39b) 



As it is clear in order to have real solutions one has to satisfy one of the following sets of conditions for the nonlinear 
coefficients: either ^22 > 912, 921 > 9iii 9ii922 > 912921 or 922 < ^12, 921 < 9iii 9ii922 < 912921- 
Then, the dark-bright type solution of Eq. ([2j) is given by 



ui{x) = (1 + acos(cjx))^/^[/i(X(x)), 
U2{x) = (l + acos(cjx))^/2/72(X(x)), 



(40a) 
(40b) 



where X{x) is given by Eq. (|2T]) . 



B. Stability of dark-bright soliton solutions 



Dark-bright solutions given by Eqs. (|4Q|) are remarkably robust to perturbations as one would in principle expect. 
All of our simulations of perturbed dark-bright solitons have lead to stable behavior. 

A couple of examples are shown in Figs. [T6l and [T71 for two different values of the nonlinearity gn = 0.26 (Fig. [T6|) 
and gii = 0.49 (Fig. [Hj). 



VI. CONCLUSIONS 



In conclusion, we have used Lie symmetries and canonical transformations to construct explicit solutions of coupled 
nonlinear Schrodinger systems with spatially inhomogeneous nonlinearities starting from those of spatially homoge- 
neous coupled nonlinear Schrodinger equations. 

Although we have restricted our attention to a few selected examples of physical relevance, the range of nonlinearities 
and potentials for which this can be done is very wide. We have constructed explicit solutions for y = with localized 
and periodic nonlinearities and discussed the parameter regimes for which solutions are stable. A similar methodology 
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FIG. 16: Dynamics of a dark-bright soliton solution given by Eq. ()40p under a perturbation of the form ([23]). Shown are (a) 
the initial density profiles of the dark (black solid) and bright (red dashed) components, and ensity plots of the time evolution 
of (b) \ui{x^t)\^ (c) |it2(x,t)|^. Parameter values are gn = 0.26, ^22 = 1, g = 0.5, a = 0.3, uu = 1, e = 0.02, q — 0.2. 
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FIG. 17: Same as in Fig. [I6]but for gxx = 0.49. 



has been applied to construct and study the stability of specific families of solutions when the potential is quadratic 
on the space coordinates. Finally we have constructed dark-bright soliton solutions that are not direct extensions of 
known solutions of the scalar spatially inhomogeneous case. 
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Appendix A: Linear stability analysis 



Let us look for perturbed solutions of Eq. ([T]) 



(Al) 
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where uqj are real solutions of the stationary problem (j2j) and \Aj \ <C uoj are perturbations. Substituting ()Aip in Eq. 
([2]) and linearizing it with respect to Aj{x) we obtain 

iAj = (A, + V{x)) Aj - V^Aj + {gjj\uojf ^ g\uoi3-j)\^) + 

9jj Wojl^Aj + 2guojUo^3-j) {^s-j + ^-j) (A2) 
Decomposing Aj into its real and imaginary parts, Aj{x^t) = Pj{x^t) + iQjix^t)^ defining the operators 

L^^^ = A^- - + '^Qjj {uojf + ^ (iio(3-j))^ ' (A3a) 
4"^ = A, - V' + gjj {uojf + ^ K(3-,))' . (A3b) 



and introducing the notation S{x) = 2^^ioi(^)'^02(^), we obtain 

-(-) 



(A4) 



Looking for solutions in the form Pj{x^t) = pj{x)e^^ and Qjix^t) = qj{x)e^^^ and defining ^7 = — k;^, Eqs. ()A4p 
reduce to 

Thus, a solitary wave solution uqj will be linearly unstable when either = Re(r2) < or Qi = lm{Q) ^ 0. 
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